
clc,clear;
close all;
warning off;
%导入数据
data =xlsread("中国主要城市坐标数据集.csv");
data=data;
% 城市坐标
cities=data;
% 城市数量
nCities = length(cities);
maxIter = 100;          % 最大迭代次数
[bP,bl,l]=ant_colony_optimization(nCities,cities,maxIter);
disp('Best path:');
disp(bP);
disp('Best length:');
disp(bl);
tusi(l,nCities,cities,bP,bl,maxIter);
%% 作图

function ktu=tusi(l,nCities,cities,bP,bl,maxIter)
l=[l(2:end)];
n=length(l);
count=1:n;
x=zeros(nCities+1,1);
y=zeros(nCities+1,1);
for i=1:nCities
    x(i)=[cities(bP(i),1)];
    y(i)=[cities(bP(i),2)];
end
x(i+1)=cities(bP(1),1);
y(i+1)=cities(bP(1),2);

figure
plot(count,l);
xlabel('迭代次数');
ylabel('最优值');
title('模拟迭代图');
hold on
figure;
plot(x,y,'-ro', 'LineWidth', 2, 'MarkerEdgeColor', 'k', 'MarkerFaceColor', 'g', 'MarkerSize', 10);
xlabel('X 坐标')
ylabel('Y 坐标')
title(['Iteration', num2str(maxIter), ' - Best Length: ', num2str(bl)])
grid on; % 打开网格
end

%% ACO(蚁群算法)
function [bestPath,bestLength,LengthT]= ant_colony_optimization(nCities,cities,maxIter)

% 距离矩阵
dist = zeros(nCities, nCities);
for i = 1:nCities
    for j = 1:nCities
        if i ~= j
            dist(i, j) = norm(cities(i, :) - cities(j, :));
        else
            dist(i, j) = inf; % 自己到自己的距离设为无穷大
        end
    end
end

% 参数设置
nAnts = 30;             % 蚂蚁数量
alpha = 1;              % 信息素重要性
beta = 5;               % 距离重要性
rho = 0.5;              % 信息素挥发系数
Q = 50;                % 信息素增加常数
alpha0=2;            %莱维飞行幂律分布的指数(0~2)
% 初始化信息素矩阵
pheromone = ones(nCities, nCities);
% 记录最佳路径及其长度
bestPath = [];
bestLength = inf;
LengthT=[];
% figure;
% 迭代优化
for iter = 1:maxIter
    % 每只蚂蚁选择路径
    paths = zeros(nAnts, nCities);
    lengths = zeros(nAnts, 1);
    for k = 1:nAnts
        paths(k, :) = construct_solution(pheromone, dist, alpha, beta, nCities);
        lengths(k) = calculate_path_length(paths(k, :), dist);
    end

    % 引入莱维飞行扰动
    for k=1:nAnts
        steps=levyfilght(alpha0);
        newpaths=prpaths(paths(k,:),steps);
        if calculate_path_length(newpaths, dist)<lengths(k)
            paths(k,:)=newpaths;
            lengths(k)=calculate_path_length(paths(k, :), dist);
        end
    end

    % 更新最佳路径
    [minLength, idx] = min(lengths);
    LengthT=[LengthT,bestLength];%作图的
    if minLength < bestLength
        bestLength = minLength;
        bestPath = paths(idx, :);
    end

    % 更新信息素
    pheromone = (1 - rho) * pheromone;
    for k = 1:nAnts
        for i = 1:nCities-1
            pheromone(paths(k, i), paths(k, i+1)) = pheromone(paths(k, i), paths(k, i+1)) + Q / lengths(k);
        end
        pheromone(paths(k, nCities), paths(k, 1)) = pheromone(paths(k, nCities), paths(k, 1)) + Q / lengths(k);
    end
    best_path=[bestPath,bestPath(1)];
    % 绘制当前最优路径
%    plot(cities(best_path, 1), cities(best_path, 2), 'o-');
%    title(['Iteration ', num2str(iter), ' - Best Length: ', num2str(bestLength)]);
%    pause(0.001); % 暂停一段时间以便观察图形变化
end
    bestPath=[bestPath,bestPath(1)];%回到起点
end
%% 选择路径
function path = construct_solution(pheromone, dist, alpha, beta, nCities)
    path = zeros(1, nCities);
    visited = false(1, nCities);
    path(1) = 5; % 选择起始城市
    visited(path(1)) = true;
    for i = 2:nCities
        currentCity = path(i-1);
        probabilities = calculate_probabilities(currentCity, visited, pheromone, dist, alpha, beta, nCities);
        path(i) = roulette_wheel_selection(probabilities);
        visited(path(i)) = true;
    end
end
%%   信息素计算
function probabilities = calculate_probabilities(currentCity, visited, pheromone, dist, alpha, beta, nCities)
    probabilities = zeros(1, nCities);
    for j = 1:nCities
        if ~visited(j)
            probabilities(j) = (pheromone(currentCity, j)^alpha) * ((1 / dist(currentCity, j))^beta);
        end
    end
    probabilities = probabilities / sum(probabilities);
end
%%  选择
function selected = roulette_wheel_selection(probabilities)
    cumulativeSum = cumsum(probabilities);
    r = rand();
    selected = find(cumulativeSum >= r, 1);
end

%% 莱维飞行步长公式
function step=levyfilght(deta)
    sigma = (gamma(1 + deta) * sin(pi * deta / 2) /(gamma((1 + deta) / 2) * deta * 2^((deta - 1) / 2)))^(1 / deta);
    u = normrnd(0, sigma);
    v = normrnd(0, 1);
    step=u/abs(v)^(1/deta);
end

%% 新路径产生
function newpath=prpaths(path,step)
    n=length(path);
    start=path(1);
    path(1)=[];%不改变起点
    newpath=path;
    % 随机选择起点
    startIdx = randi(n-1); 
    % 生成莱维步长        
    stepLength = ceil(step);
    % 确定交换的范围
    endIdx = mod(startIdx + stepLength - 1, n-1) + 1;
    % 交换路径中的两个城市
    temp = newpath(startIdx);
    newpath(startIdx) = newpath(endIdx);
    newpath(endIdx) = temp;
    newpath=[start,newpath];
end
%%  目标值计算
function length0 = calculate_path_length(path, dist)
    length0 = 0;
    nCities = length(path);
    for i = 1:nCities-1
        length0 = length0 + dist(path(i), path(i+1));
    end
    length0 = length0 + dist(path(nCities), path(1)); % 回到起点
end

